First, some clarification: in the context of, typically, heat
propagation and/or (magneto)hydrodynamics along a line one pays
attention to the difference boundary-value problem  on a lattice
or ``grid" or ``net".

The recommended trick is that one reduces the second-order
equation for intensities $y(x)$ on the grid to the system of the
first-order equations for $y(x)$ and the first differences or
``flows" $w(x)$. At generalized boundary conditions one rotates
the intensities $y(x)$ and flows $w(x)$ into certain
``orthogonally rotated" two-dimensional vectors $\vec{s}(x)$ with
evolution controlled by the matrices which are lower triangular
for certain grid-point-dependent rotation angles $\alpha(x)$.

One ends up with the algorithm which solves the problem as
decoupled recurrences ``from left to right" for the upper
components of $\vec{s}(x)$ and then, similarly but in the opposite
direction, for the lower components of $\vec{s}(x)$. Recommended
by the authors whenever the ``flows" play a central role since
their solution method incorporates them directly. The use of
rotations may also make the stability achievable under reasonably
weakened constraints.


MR2164578
Monastyrny\u\i, P. I. ; Shvakel\cprime, A. I. A data-driven
orthogonal sweep method for difference equations with strongly
varying coefficients. (Russian) Vests\=\i Nats. Akad. Navuk
Belarus\=\i Ser. F\=\i z.-Mat. Navuk 2004, no. 3, 22--28, 126.